home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / rwvector.lha / RWVector2.1 / src / dfft.cc < prev    next >
C/C++ Source or Header  |  1989-08-18  |  5KB  |  167 lines
  1. /*
  2.  *    Definitions for the double precision FFT server.
  3.  *
  4.  *    Copyright (C) 1988, 1989.
  5.  *
  6.  *    Dr. Thomas Keffer
  7.  *    Rogue Wave Associates
  8.  *    P.O. Box 85341
  9.  *    Seattle WA 98145-1341
  10.  *
  11.  *    Permission to use, copy, modify, and distribute this
  12.  *    software and its documentation for any purpose and
  13.  *    without fee is hereby granted, provided that the
  14.  *    above copyright notice appear in all copies and that
  15.  *    both that copyright notice and this permission notice
  16.  *    appear in supporting documentation.
  17.  *    
  18.  *    This software is provided "as is" without any
  19.  *    expressed or implied warranty.
  20.  *
  21.  *
  22.  *    @(#)dfft.cc    2.1    8/18/89
  23.  */
  24.  
  25.  
  26. #include "rw/DoubleFFT.h"
  27. #include "rw/mathpack.h"
  28.  
  29. DoubleFFTServer::DoubleFFTServer()
  30. {
  31.   server_N=0;
  32. }
  33.  
  34. DoubleFFTServer::DoubleFFTServer(unsigned N)
  35. {
  36.   setOrder(N);
  37. }
  38.  
  39. DoubleFFTServer::DoubleFFTServer(const DoubleFFTServer& s)
  40.      : (s), roots_of_1(s.roots_of_1), conjroots_of_1(s.conjroots_of_1)
  41. {
  42.   server_N = s.server_N;
  43. }
  44.  
  45. void
  46. DoubleFFTServer::operator=(const DoubleFFTServer& s)
  47. {
  48.   DComplexFFTServer::operator=(s);
  49.   server_N = s.server_N;
  50.   roots_of_1.reference(s.roots_of_1);
  51.   conjroots_of_1.reference(s.conjroots_of_1);
  52. }
  53.  
  54. void
  55. DoubleFFTServer::setOrder(unsigned N)
  56. {
  57.   server_N = N;
  58.   DComplexFFTServer::setOrder(N);
  59.   roots_of_1.reference(rootsOfOne(2*N, N/2+1));
  60.   conjroots_of_1.reference(conj(roots_of_1));
  61. }
  62.  
  63. /*
  64. Performs DFT of a real sequence.  Must have an even number of points.
  65. The forward fourier transform of a real sequence is a complex
  66. conjugate even sequence.  If the original sequence had 2N points, the
  67. result will have N + 1 complex points, for a total of 2N+2 points.
  68. The extra two points are the imaginary parts of C(0) and C(N).  Both
  69. are always zero.  Reference: Cooley, et al.  (1970) J. Sound Vib.
  70. (12) 315--337.  This is algorithm #4.  See the class description
  71. header for more information. 
  72. */
  73.  
  74. DComplexVec
  75. DoubleFFTServer::fourier(const DoubleVec& v)
  76. {
  77.   unsigned Nstore    = v.length();
  78.   unsigned N        = Nstore/2;
  79.   unsigned Nhalf    = N/2;
  80.   unsigned Nhalfp1    = N/2+1;
  81.   fortran_int Nint    = N;
  82.  
  83.   // Various things that could go wrong:
  84.   checkEven(Nstore);
  85.   if(!Nstore) return DComplexVec();
  86.  
  87.   if(N != server_N) setOrder(N);
  88.  
  89. #if COMPLEX_PACKS
  90.   // Bit of a cludge, in the interests of efficiency:
  91.   DComplexVec A((DComplex*)v.data(), N);
  92.   DCfftf_(&Nint, A.data(), the_weights.data());
  93. #else
  94.   DComplexVec tempcomplex(v.slice(0,N,2), v.slice(1,N,2));
  95.   DComplexVec A = DComplexFFTServer::fourier(tempcomplex);
  96. #endif
  97.  
  98.   DComplexVec results(N+1);
  99.   DComplexVec A1(Nhalfp1);
  100.   DComplexVec A2(Nhalfp1);
  101.  
  102.   A1.slice(1,Nhalf,1) = conj(A.slice(N-1,Nhalf,-1)) + A.slice(1,Nhalf,1);
  103.   A2.slice(1,Nhalf,1) = DComplexVec(imag(A.slice(N-1,Nhalf,-1))+imag(A.slice(1,Nhalf,1)),
  104.                     real(A.slice(N-1,Nhalf,-1))-real(A.slice(1,Nhalf,1)));
  105.  
  106.   A1(0) = DComplex(2.0*real(A(0)),0);
  107.   A2(0) = DComplex(2.0*imag(A(0)),0);
  108.  
  109.   DComplexVec temp = A2 * conjroots_of_1;
  110.   results.slice(0,Nhalfp1, 1) = DComplex(0.5) *     (A1 + temp);
  111.   results.slice(N,Nhalfp1,-1) = DComplex(0.5) * conj(A1 - temp);
  112.  
  113.   return results;
  114. }
  115.  
  116. /*
  117. Take the IDFT of a complex conjugate even sequence.  V is assumed to
  118. be the lower half of the complex conjugate even sequence.  See the
  119. comments in the class header about the format of even sequences.  The
  120. results is a real periodic sequence.  Reference: Cooley, et al.
  121. (1970) J. Sound Vib.  (12) 315--337.  This is algorithm #5.
  122. */
  123.  
  124. DoubleVec
  125. DoubleFFTServer::ifourier(const DComplexVec& v)
  126. {
  127.   unsigned N        = v.length()-1;
  128.   unsigned Ntot        = 2*N;
  129.   unsigned Nhalf    = N/2;
  130.   unsigned Nhalfp1    = N/2+1;
  131.  
  132.   if(N != server_N) setOrder(N);
  133.  
  134.   DComplexVec lowslice = v.slice(0,Nhalfp1,1);
  135.   DComplexVec hislice = conj(v.slice(N,Nhalfp1,-1));
  136.  
  137.   DComplexVec A1 = lowslice + hislice;
  138.   DComplexVec A2 = (lowslice - hislice) * roots_of_1;
  139.   
  140.   DComplexVec A(N);
  141.   A.slice(0,Nhalfp1,1)  = DComplexVec(real(A1)-imag(A2), real(A2)+imag(A1));
  142.   A.slice(N-1,Nhalf,-1) = DComplexVec(real(A1.slice(1,Nhalf,1))+imag(A2.slice(1,Nhalf,1)),
  143.                       real(A2.slice(1,Nhalf,1))-imag(A1.slice(1,Nhalf,1)));
  144.   
  145.   DComplexVec X = DComplexFFTServer::ifourier(A);
  146.  
  147. #if COMPLEX_PACKS
  148.   DoubleVec results((double*)X.data(), Ntot);
  149. #else
  150.   DoubleVec results(Ntot);
  151.   results.slice(0,N,2) = real(X);
  152.   results.slice(1,N,2) = imag(X);
  153. #endif
  154.   return results;
  155. }  
  156.  
  157. void
  158. DoubleFFTServer::checkEven(int l)
  159. {
  160.   if( l%2 ){
  161.     char msg[60];
  162.     sprintf(msg, "Length (%d) of a real series must be even.", l);
  163.     RWnote("DoubleFFTServer::[i]fourier()", msg);
  164.     RWerror(DEFAULT);
  165.   }
  166. }
  167.